Viscous potential flow analysis of peripheral heavy ion collisions 
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The conditions for the development of a Kelvin-Helmholtz Instability (KHI) for the Quark-gluon 
Plasma (QGP) flow in a peripheral heavy-ion collision is investigated. The projectile and target 
side particles are separated by an energetically motivated hypothetical surface, characterized with 
a phenomenological surface tension. In such a view, a classical potential flow approximation is 
considered and the onset of the KHI is studied. The growth rate of the instability is computed as 
function of phenomenological parameters characteristic for the QGP fluid: viscosity, surface tension 
and flow layer thickness. 
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I. INTRODUCTION 

The first models of high energy heavy ion collisions in 
the 1970s were successful by assuming highly idealized 
shock fronts where the matter was heated up and com- 
pressed in a front (having a discontinuity in perfect fluid 
flow [IlHj)- This led to high pressured, shock compressed 
domains which collectively deflected the incoming nuclear 
fluid. The observation of this directed flow (side splash 
of bounce-off) was the first proof of the collective fluid 
dynamical behaviour of nuclear matter [3J. Recent the- 
oretical developments and experimental observation of 
high multipolarity fluctuations indicate that the QGP is 
a low viscosity fluid, which makes turbulent phenomena 
possible 0HB] 
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FIG. 1. Sketch of a collision: (a) is a view in the transverse, 
[x,y] plane and (b) is an illustration in the reaction, [x,z] plane. 
The almond shape in the middle of figure (a) is the participant 
zone of the event. Right after the collision, streaks are formed 
and the top streaks move along the z direction while bottom 
ones move along the — z direction. Due to this velocity shear, 
an instability wave will appear on the interface plane between 
the top and bottom sheets. 



Here we adopt again a fluid dynamical picture and dis- 
cuss the strong shear flow arising in the initial states of 
peripheral heavy ion collisions at ultra-relativistic ener- 
gies, which may lead to KHI under favorable conditions, 
as discussed recently [4] . 

A simple analytic study showing the development of 
the KHI in a highly idealized situation is discussed. The 



investigated phenomenon has some resemblance to the 
initial state of a peripheral heavy ion collision. In these 
collisions the collective flow should be a " shear flow" be- 
cause the top participant layers, move nearly with pro- 
jectile velocity while the bottom layers with the target 
velocity. 

In the reaction plane, the height of the participant pro- 
file is L = 2R — b, where b is the impact parameter, and 
the half height is I = L/2 = (2R - b)/2. In the follow- 
ing we will denote by n the nuclear matter density and 
by r] its phenomenological viscosity. The matter coming 
from different sources is marked by subscript/superscript 
t (top) and b (bottom), see Fig. [I] 

For simplicity reasons we consider a two dimensional 
non-relativistic dynamics in the reaction plane. The po- 
sition vector is x = (x,z), and the velocity vector is 
v = (v x ,v z ). 
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FIG. 2. The velocity along the z axis, is represented by the 
dotted curve, calculated in our CFD model and presented in 
Ref. The full line with two singularities in its derivative 
is the case used in Ref. [7]. Here we have to mention that 
the velocity profile of the dotted curve will induce the KHI 
effect, while the velocity profile illustrated with the dashed 
curve will not, see chapter 8 of Ref. [8]. 

The nearly perfect QGP provides a possibility of a 
strong idealization in this situation. The velocity profiles 
for v z presented on Fig. [2] illustrates (the dotted line) 
that the KHI develops if we have a strong shear flow at 
the x — plane, which leads to large vorticity and circu- 
lation. With decreasing viscosity we could idealize this 
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configuration in a way that the vorticity is constrained 
into a narrow layer around the x = plane, while the 
circulation remains constant. In a limiting case if we 
constrain the vorticity and shear to the dividing, x = 0, 
plane, this plane would represent an infinite vorticity, 
providing the same circulation for a trajectory surround- 
ing the dividing plane. Q In this limiting case the velocity 
profile is idealized to the one indicated by the full line in 
Fig. [2] Thus in the top and bottom domain the flow has 
no shear, and can be described as potential flow, which is 
an important simplification and idealization. This makes 
the analytic study of the KHI possible. 

II. THE ANALYTIC MODEL 

Thus, the idealized dividing layer represents a disconti- 
nuity of the flow velocity (i.e. unconstrained slip). At the 
same time a small viscosity would contribute the trans- 
verse momentum transfer to a small distance across the 
dividing front. The particles in this narrow layer, scatter- 
ing over from the other side of the dividing plane, would 
have a high relative velocity, so this layer would exhibit 
an extra energy increase compared to the general fluid 
body on the top or the bottom side of our system. This 
can be taken into account as an effective surface energy of 
the dividing layer. This surface energy can be estimated 
both from a microscopic kinetic theory approach, or from 
a rough energy balance calculation. In a microscopic ap- 
proach one would assume that the this extra energy de- 
pends on the (viscosity dependent) thickness of the layer 
and the (temperature dependent) rate of transverse flow 
crossing the dividing plane. A quantitative estimate of 
this surface energy in this idealized situation is not fea- 
sible, but its existence and a rather qualitative estimate 
can be made. In contrast with this, a phenomenological 
energy balance calculation is easier to perform. The ad- 
vantage of such an approach would be that one does not 
rely on further estimates for the involved physical param- 
eters. Here, we will use this later method to approximate 
the surface tension of the dividing layer. 

Following ref. [7] we idealize the problem and assume 
an initial state where the shear is localized at the dividing 
plane between the top (t) half and the bottom (b) half of 
the fluids, in order that we can use the potential flow de- 
scription in the top and bottom parts of the fluid, see Fig. 
[3j We assume that the fluid in the top and bottom parts 
are allowed to slip at the top and bottom boundaries as 
well as at the dividing surface between them. We will ref- 
erence these as unconstrained slip-conditions. The initial 
flow velocity is assumed to be uniform in the two layers, 
so that for the top layer t>* = (0, Ut) for < x < I and 



1 The conservation of circulation occurs in classical, barotropic 
flow. In QGP the temperature dominates the pressure change, 
so the circulation is not conserved but decreases during the ex- 
pansion of the system [9]. 



for the bottom layer v b = (0, Ub) for — I < x < initially. 
This means that initially the amplitude of the wave-like 
instability is extremely small, and we are looking for the 
conditions to have a growing amplitude for this instabil- 
ity. For the sake of analytic model we assume that the 
density is constant. Numerical studies |4j show that this 
constraint can be relaxed. 
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FIG. 3. The profile of the top and bottom fluids layers with a 
dividing surface wave on it. The top external fluid moves with 
velocity Ut along the z-direction, while the bottom fluid moves 
with velocity Ub along —z-direction. The slip is unconstrained 
at the dividing surface. 



Just as in ref. [3] we will use an average energy and 
mass density for our estimates. Since the effective parti- 
cle density is constant the continuity equation of the flow 
velocity, v, will become: 

V-v = 0. (1) 

For the top and bottom parts of the fluid we assume 
small velocities and neglect velocity gradients and we as- 
sume that the rotation of the flow velocity is V x v = 0. 
Under such conditions we can describe the flow as a 
potential flow, i.e., v = V0, where (f> is the velocity 
potential, cf> = (f>(t,x,z). The continuity equation is 
A(f> = \7 2 (j) = 0, applied this for the top and bottom 
layers gives: 

V 2 <j)t = for < x < I, (2) 
V 2 <j) b = for -I < x < 0. (3) 

We assume that due to the raising instability the ini- 
tially plane interface will experience a perturbation and 
will deviate from the x = plane. The height of the devi- 
ation from the x = plane is denoted by ft, = h(t, z), and 
it is taken as wave-like perturbation in the z direction 
with wave-number k. We also allow the amplitude at a 
given coordinate to change in time. The most general 
form for such a wave-like interface would be 

h(t,z) =A e( at+lkz \ (4) 

where, a is the complex growth-rate, and Aq is a complex 
amplitude. 

Consequently, the fluid on the top and bottom sides 
next to the dividing surface will have vertical velocity 
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components: 



vt; b = — = 



dh 
dt 



dh 

dt 



dh 

' dz ' 



(5) 



where Ut and Ub are the flow velocities of the fluid at 
the top and bottom of the bounding surfaces, respec- 
tively. These are assumed to be the average velocity of 
the fluids at the surface neglecting the horizontal velocity 
fluctuations arising from wave formation. Initially these 
velocity fluctuations are small, and our aim is to study 
the initial development of KHI. At the h(t, z) dividing 
surface we also assume unconstrained slip conditions of 
the two inversely flowing fluid slabs. 

The boundary conditions for the external border of the 
profile are: 



„t,b 



dx 



at x = ±L 



Initially (at time t = 0), one would also have to satisfy: 
v t /(t = 0,x = ±l) = U t , b (7) 



At each time moment the potential (f>t and <pb is satis- 
fying Eq. ( 2|3|5|6 ) for the top and bottom sides respec- 
tively: 
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(8) 



9a; 2 
dh dh 

" x ~dt = lh 
t,b d $t.b 
V * = ^x- 

Assuming for the interface, h{t,z), in Eq. Q, a wave- 
like perturbation, which is symmetric in ±x and expo- 
nentially decreasing away from the surface, the solution 
can be searched in the form: 

4>t,b = A ttb cosh[fc(z - 0]e (<T * +ifcz) + zU t , b , (9) 

where At, A b and Aq are the complex amplitudes, a is 
the growth rate and k is the wave number. From the 
kinematic conditions on the dividing layer Eq. ([5]), we 
get the following equations at the dividing surface: 



(<r + ikU t ,b)A = TkA tjb sinh(fc^) 



(10) 



The pressure (p), viscosity (77) and surface tension (7) 
balance at the interface writes as: 

„ dvt . „ dvi , d 2 h , . 

The surface energy and consequently the surface tension 
of the dividing layer will be approximated later. As it 
was already emphasized in the introductory paragraphs, 
although the top (t) and bottom (b) sides are of the same 
nuclear matter, the velocity jump or the sharp velocity 
change contribute to additional surface energy due to the 
large shear at the interface exhibiting extra energy or to a 
smaller extent by the momentum dependance of nuclear 
interaction potential. 



Since we have unconstrained slip conditions on the di- 
viding surface between the top an bottom layer, p t and 
p b can be written by the classical equation of motion 
without the viscous term as: 



dv 



t,b 



u, 



do 



t,b 



dt ' "*' b dz ' dz ' 

Then, first we apply V z on both sides of the equation and 
substitute equation of continuity, d z v z — —d x v x , into it: 



Op, 



Lb 



(12) 
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„i,6 



P( 



dtdx 



l.b 



dxd: 
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d 2 p t ,, 
dz 2 



(13) 



Here p is the effective mass density of the QGP, we use 
p = 10 GeV/fm 3 c 2 [1] in our work. In order to substi- 



tute the above equations into Eq. (11), we consider the 
second order derivative of Eq. (11) as a function of z, 



(g) and substitute Eq. (13) into it. Thus the pressure, vis- 



cosity and surface tension balance will be written in the 
following form: 

£)2„,t A2„,t Q3 v t 



dtdx 



U, 



dxdz 
d 3 vt 



dxdz 



d A h 



(14) 



^ dtdx b dxdz^ ^ dxdz 2 ^ dz A 

By inserting the velocity derived from Eq. ^ and 
the considered interface profile, Eq. into the above 
equation, and expressing the top and bottom amplitudes, 
A t b from Eq. ( 10 ), after simplifying all over with A , and 
putting the condition x = 0, we obtain an equation for a 
and k: 



[p {a + ikU t ) 2 

+ [p(<7 + ikU b ] 

+ jk 3 = 0. 



■f 2r]k 2 ((T + ikUt)} coth(fcZ) 
2 + 2rjk 2 (a + ikU b )\ coth(fcJ) 



(15) 



Considering this as an equation for a, one can write it in 
a simplified form as 



Aa 2 + 2Ba + C = 0, 
where the coefficients, A, B, C are defined as: 
A = 2p coth(fc/), 



(16) 



ikp {Ub+Ut) coth(kl) 



B = 2k 2 i] coth(fcZ) 

= B R +iB l7 
C = -k 2 p coth(kl)(U 2 + U 2 ) + 

U b ) = 



(17) 



+ 2ik 3 r] coth(AZ) (U t 
The solution is 
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D R = k 2 p 2 coth 2 (kl)(U t - U b ) 2 

+ 4r] 2 fc 4 coth 2 (kl) - 2p coth(kl)-/k 3 , 
D T = 0, 



(18) 



(19) 
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thus the real part and the imaginary part can be ex- 
pressed as: 
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(20) 
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FIG. 4. (color online) The real part of the growth rate, ctr, 
is shown as function of the viscosity rj. The full (red) line 
is for the surface tension 7 = 0.4 GeV/fm 2 and the dashed 
(blue) line is for 7 = 3.5 GeV/fm 2 . The wave number, fc, 
is taken to be k — 0.6 fm" 1 and the effective mass density 
is p — 10 Gev/fm 2 c 2 . The growth rate decreases when the 
viscosity increases suggesting that the KHI grows weaker for 
a more viscous fluid. 



In heavy ion collisions, the matter will expand after 
the collision, and in fact there is no external boundary 
(top and bottom) of the fluid shown in Fig. [3] If we 
assume I — > 00, the above equations can be simplified as: 



e>7 = 



± 



k 2 rj 

P V p " 
k(U t + U b ) 



'&V k 2 (U t -U b ) 2 7 fc 3 



2p 



(21) 



(22) 



In Eq. ([21]), for the typical parameters of a peripheral 
heavy ion collision, the real part of the growth rate, <tr, 
is dominantly dependent on the viscosity 77, namely the 
first term and the first term in the square root. In our 
expanding system the dominant wave number of KHI is 
changing with time. 



III. RESULTS 

According to the CFD observations [3], initially we 
have a small wave formation with k 1 fm -1 , but with 
time and expansion, the possible largest wave length 
takes over with k f=s 0.6 fm -1 , which decreases fur- 
ther with the expansion of the system. By assuming 
\Ut — Ub\ = 0.8 c, we can obtain the growth rate de- 
pendence of the viscosity, 77, and the surface tension, 7, 
which are shown in Figs. [4]and[5j 




7 (GeV/fm-) 

FIG. 5. (color online) The real part of the growth rate, o\r, 
as a function of the surface tension 7 with different values 
of viscosity 77. The full (red) line represents 77 = 3 GeV/fm 2 
and the dashed (blue) line represents 77 = 6 GeV/fm 2 . The 
wave number k is taken as 0.6 fm" 1 and the effective mass 
density p is 10 Gev/fm 2 c 2 . As we can see in the figure, 
the two curves cross each other at or = 0, which is around 
7 = 5.3 GeV/fm 2 , and then the growth rate becomes negative. 
With bigger surface tension the KHI effect is less probable to 
appear. 



Similarly to the viscosity the effective surface energy 
also influences the growth rate of KHI. As expected larger- 
surface tension or surface energy damps the growth of 
KHI. Beyond a critical surface energy (in our model at 
jcrit ~ 5.3 GeV/fm 2 ) the surface tension will lead to a 
decrease in the KHI. Interestingly this threshold value is 
independent of the viscosity. This is part of the general 
feature that the behavior of the zero-growth (a R = 0) 
curve is independent of the value of the viscosity in this 
model. The growth rate and damping rate are of course 
depend on the viscosity. 

The condition to have a growing instability is to have 
a solution with <jr > 0. Taking into account that D is a 
real number (Dj = 0, D R > 0), and B R > 0, from ( [20] ) 
it follows that in order to have a positive growth rate, 
(or, > 0), one has to satisfy the condition: 



(23) 



Thus, using (17) and (19) we get the condition for posi- 
tive growth: 



V > 



2jk 



p coth(fci) 



(24) 



where V = U t — Ub- 

The above condition will limit the region of the (V, k) 
parameter space where the KHI can evolve. One should 
also keep in mind the results obtained in [4], regarding 
the acceptable wave numbers, k, for the considered wave- 
like instability. Definitely there is a lower cutoff (k m i n ) 
governed by the beam-directed longitudinal length of the 
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flow, L 



•2tt 



(25) 



For the b = 0-5b max and b = 0.7b max impact parameter 
values the calculations in [3] leads to k min = 0.598 fm -1 
and kmin = 0.479 fm -1 values, respectively. There is also 
an upper limit for the wave-numbers, k max governed by 
the Kolmogorov length scale, Xk ■ 



h — 



2tt 
X~ 



K 



(26) 



According to [4] this characteristic length-scale is esti- 
mated for the above given impact parameters as: Xk ~ 
3.5 fm and Xk ~ 2.5 fm, leading to k max = 1.79 fm -1 
and k m ax = 2.51 fm -1 values, respectively. 

For the peripheral Pb+Pb collisions, the radius of Pb is 
R = 7 fm, thus bmax — 14 fm. In order to get the param- 
eter space where the KHI will growth let us estimate now 
the value of the surface tension. As it has been discussed 
in the introductory part this surface energy comes from 
the energy excess of the unbalanced energy flow in the 
two layers. Although a theory based on kinetic considera- 
tions would capture more from the involved physics, here 
we just consider a simple approach based on the energy 
balance. The reason for doing this is that less number of 
phcnomenological parameters are needed. 
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(a) 



(b) 



FIG. 6. Two different velocity profiles: (a) is the profile used 
in our present work, it has two distinct layers with two con- 
stant velocity Ut and Ub, while (b) has a flow transition from 
the two layers, and at the diving surface the velocity is small- 
est. 

The flow assumed in the present work has a perpen- 
dicular velocity profile illustrated in Fig. [6^,. This means 
that a smooth velocity profile (Fig. [6Jd), characterizing 
a stable and balanced viscous flow is not formed. In the 
case illustrated in Fig. [6^i one would assume that there 
are two distinct layers flowing with velocities U t and U - 
For the balanced flow illustrated in Fig. [SJd, one would 
observe a smooth flow velocity transition from Ut to U - 
It is obvious that in the laboratory frame, this later flow 
has less kinetic energy in the z direction than the previ- 
ous one. The difference between the two kinetic energies 
can be accounted as the energy surplus of the dividing 
layer. If we denote the contact surface between the flows 
in the top and bottom layers by S, the surface tension 
could be estimated as 



E L 



E kz 



(27) 



where E^' z b denotes the kinetic energy of the flow in the 
z direction for the profile illustrated in Fig. [6^, and|6]3, 
respectively. The total relativistic kinetic energy of the 
system, in the laboratory frame is 



E k = 2M Pb c 




(28) 



where V = Ut — Ub is the relative speed of the two pro- 
jectiles, and Mph is the mass of the collided Pb ions. 
Assuming that the participating zone in the collision has 
a surface q x irR 2 , (the overlapping regions are only a q 
part of the possible ones) and the kinetic energy of the 
participating particles in this zone is distributed equally 
in all the directions of the space, a rough approxima- 
tions for E^ z would be: E% z = q_j^- On the other hand, 
for the flow illustrated in Fig. [6}d, due to the balanced 
velocity profile, a part of this kinetic energy has to be dis- 
sipated, and assuming a linear velocity profile one gets: 
E h kz = 1/2 E% z . The above arguments lead us to a first 
approximation of the surface tension value: 



q MpbC 




(29) 



Assuming q sa 0.5 and estimating the surface of the 
dividing layer, S, from [3], one gets the values of 7 for 
different impact parameter values. 

The surface tension is estimated to be 7 = 
0.4 GeV/fm 2 from Eq. (29). This value is used in the 
following examples. The critical velocity Eq. ( p4| for 
different impact parameters is shown in Fig. [7] These 
curves show the border of instability of the growth rate, 
up. The curves divide the space into two areas, the upper 
side above the curve is the region where the instability 
grows and the area below the critical velocity curve is 
where the instability does not grow. The KHI develop- 
ment region is also limited by the k m m and k rnax values 
as drawn in figure Fig. [7] 

The above consideration is for or — 0, however, this 
does not show the ^-dependence of the growth. In order 
to see how the instability depends on the viscosity, ?y, we 
can cast Eq. ( 20 1 into the form: 



k 2, q 
P 




coth(M) J 



(30) 



This suggests that with our characteristic parameters the 
dependence on the thickness of the fluid layer, I, is weak 
as shown in Fig. [8] 
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FIG. 7. (color online) The critical velocity V = Ut — Ub 
as function of the wave number, at the condition of van- 
ishing growth rate, <jr — 0. The red full line, blue dashed 
line and black dot-dashed line are for impact parameters 
b = 0.5,0.7,0.9 b m ax, respectively. On the graph we also 
illustrated the two natural boundaries fc min and k max for 
b — 0.5 b m ax- The KHI will evolve thus above the critical 
velocity curves and between these two limits. For increasing 
impact parameters, the instability is less able to grow and the 
system tends to be stable. 
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FIG. 8. (color online) The growth rate or as a function of the 
viscosity rj at different values of I, black line is for I = 0.1 fm, 
the red dashed line is for I — 2 fm and the blue dot-dashed 
line is for I = oo. The wave number is k — 0.6 fm -1 , the 
surface tension is 7 = 0.4 Gev/fm 2 , the relative velocity is 
V = 0.8 c and p = 10 Gev/fm 2 c 2 . The growth rate depends 
weakly on while it depends significantly on the viscosity, 
increasing strongly for small viscosity values. 

IV. CONCLUSIONS 

In classical gravitational water waves the wave forma- 
tion and wave speed depends strongly on the depth of the 



water, i.e. the layer thickness, I. In a heavy ion collisions 
the role of the layer thickness is different. The material 
properties of the top and bottom layers are not different, 
these are separated from another by the relatively thin 
layer of large shear. Still the occurrence of KHI in such 
conditions is not uncommon as it is frequently observed, 
as turbulence during airplane flights, or it is even visible 
if the air has high humidity and the condensation makes 
the KHI visible. 

In peripheral heavy ion collisions the layer thickness is 
given in the initial state, but there is no solid boundary 
and the system expands in all directions. Thus, for this 
physical situation the large or infinite layer thickness is 
more relevant in this model, even if the initial layer thick- 
ness is finite and usually smaller than the longitudinal 
size of the initial state. 

Large viscosity or the corresponding low Reynolds 
number prevent the development of turbulence and KHI, 
so that these phenomena appear only above a critical 
Reynolds number. This critical Reynolds number de- 
pends on the flow configuration, so it is separately anal- 
ysed for the KHI also, see ref. |4j. The present study 
confirms that the dependence of the growth rate on the 
viscosity reflects the usual tendency that instability and 
turbulence increases with smaller viscosity. 

When the KHI develops between two fluids (e.g. 
air/water or air/oil) the large surface tension difference 
at the interface damps the development of the instability, 
this is well known for sailors for centuries. If KHI devel- 
ops inside one fluid, like in air or in quark gluon fluid, 
there is no surface tension in the classical sense, but the 
layer with large shear has extra energy, and it leads to an 
effective surface tension, which hinders the development 
of KHI. 

We presented a strongly idealized analytic model for 
the development of the Kelvin-Hclmholtz Instability in 
ultra-relativistic heavy ion reactions. We compressed the 
shear zone into a central infinitesimal layer, following the 
idea of ref. [7], and assumed that the remaining flow 
can be approximated as potential flow. The idealized di- 
viding layer was attributed a surface energy and uncon- 
strained slip between the top and bottom fluid layers. It 
is interesting that in this model the KHI is developing 
under similar conditions, as in numerical high resolution 
relativistic fluid dynamical calculations 0]. This model 
also shows that critical size KHI may occur for low vis- 
cosity QGP. 
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